Orbital magnetization and Chern number in a supercell framework: 

Single k-point formula 
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The key formula for computing the orbital magnetization of a crystalline system has been recently 
found [D. Ceresoli, T. Thonhauser, D. Vanderbilt, R. Resta, Phys. Rev. B 74, 024408 (2006)]: it 
is given in terms of a Brillouin-zone integral, which is discretized on a reciprocal-space mesh for 
numerical implementation. We find here the single k-point limit, useful for large enough supercells, 
and particularly in the framework of Car-Parrinello simulations for noncrystalline systems. We 
validate our formula on the test case of a crystalline system, where the supercell is chosen as a large 
multiple of the elementary cell. We also show that — somewhat counterintuitively — even the Chern 
number (in 2d) can be evaluated using a single Hamiltonian diagonalization. 

PACS numbers: 75.10.-b, 75.10.Lp, 73.20.At, 73.43.-f 
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The position operator r is ill-defined within peri- 
odic boundary conditions. Owing to this, both the 
macroscopic (electrical) polarization and the macro- 
scopic orbital magnetization are nontrivial quantities in 
condensed-matter theory. The former has been suc- 
cessfully tamed since the early 1990s, when the mod- 
ern theory of polarization, based on a Berry phase, was 
developed. 1,2 The latter, instead, remained an unsolved 
problem until 2005. Since then, several important papers 
have appeared 3-7 and continue to appear. 8 Before 2005 
only linear-response properties related to orbital magne- 
tization were successfully addressed, 9,10 while we stress 
that the present work, as well as Refs. 3-8, addresses 
"magnetization itself" . 

A general formula, valid for crystalline systems within 
a given single-particle Hamiltonian, was provided in 
Ref. 6, hereafter referred to as I. This is the magnetic 
analogue of the (by now famous) King-Smith and Van- 
derbilt formula for electrical polarization. 1 Both formu- 
las are discretized on a regular mesh of k points for 
numerical implementation. However, most simulations 
for noncrystalline systems, particularly those of the Car- 
Parrinello type, 11 are routinely performed by diagonaliz- 
ing the Hamiltonian at a single k point (the T point) in a 
large supercell. The reduction from many points to a sin- 
gle point is far from being trivial; nonetheless a successful 
single-point formula for electrical polarization emerged 
since 1996, and is universally used since then. 12,13 We 
provide here the magnetic analogue of such formula. As 
a byproduct, one can even evaluate the Chern number 14 
using a single k point. This looks like an oxymoron, given 
that the Chern number is by definition a loop integral in 
reciprocal space: but our formula can be regarded as the 
limiting case where the loop shrinks to a single point. 

The general formula of I applies to normal periodic 
insulators (where the Chern invariant vanishes), Chern 
insulators (where the Chern invariant is nonzero), and 
metals. Aiming at first-principle implementations within 



any flavor of DFT, the single-particle Hamiltonian is the 
Kohn-Sham one. 15 As for the analogous case of electrical 
polarization, there is no guarantee that the Kohn-Sham 
magnetization coincides with the actual one. Nonethe- 
less, previous studies within linear-response methods in- 
dicate that even for orbital magnetization the error is 
small. 10 

As in I, we assume a vanishing macroscopic magnetic 
field B, hence a lattice-periodical Hamiltonian. We let 
e n k and |^ n k) be the Bloch eigenvalues and eigenvectors 
of H, respectively, and w„k(r) = e~ ik ' T ipnk( r ) be the cor- 
responding eigenfunctions of the effective Hamiltonian 



H k = c 



k '7/e' k r : 



(1) 



we normalize them to one over the crystal cell of volume 
V. As in I, the notation is intended to be flexible as 
regards the spin character of the electrons. If we deal 
with spinless electrons, then n is a simple index labeling 
the occupied Bloch states; factors of two may trivially 
be inserted if one has in mind degenerate, independent 
spin channels. For the sake of simplicity, we rule out 
the metallic case here. For both normal insulators and 
Chern insulators the macroscopic orbital magnetization 
M is, according to I: 



M 7 = 



2c(2tt) 3 



(2) 
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where Greek subscripts are Cartesian indices, £ ia p is the 
antisymmetric tensor, d a = d/dk a , /i is the chemical po- 
tential (Fermi energy), the integration is over the recip- 
rocal cell, and the sum over Cartesian indices is implied. 
For insulators the number of states with energy e„k < M 
is independent of k; we implicitly understand the sum in 
Eq. (2) as on these states only. 
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As usual, a noncrystalline system can be dealt with in 
a supercell framework, by addressing an artificial crystal 
with a large unit cell, actually larger than the relevant 
correlation length. One key feature of Eq. (2) is its gauge- 
invariance in the generalized sense: by this we mean that 
Eq. (2) is invariant by unitary mixing of the occupied 
states among themselves at a any given k. Thanks to 
such key feature Eq. (2) is invariant by cell doubling. In 
fact, starting with a cell (or supercell) of given size, we 
may regard the same physical system as having double 
periodicity (in all directions), in which case the integra- 
tion domain in Eq. (2) gets "folded" and shrinks by a 
factor 1/8, while the number of occupied eigenstates gets 
multiplied by a factor of 8. It is easy to realize that 
these are in fact the same eigenstates as in the unfolded 
case, apart possibly for a unitary transformation, irrel- 
evant here. As for actual computations, the discretized 
form of Eq. (2) adopted in I turns to be invariant by cell 
doubling to within 10~ 5 , provided the k-point mesh is 
chosen consistently. 

For a large supercell of volume V the integration do- 
main in Eq. (2) becomes small. Therefore the integral 
can be accurately approximated by the value of the inte- 
grand at k = 0, times the reciprocal-cell volume: 

M 7 ~ -7^7 Im ^2(d a u n0 \(H + e„ - 2yu) \dpu n0 ) . 

n 

(3) 

Notice that Eq. (2) can be safely approximated with 
Eq. (3) because the integrand is a gauge-invariant quan- 
tity; the apparently analogous case of polarization is 
more difficult, since the integrand therein is gauge- 
dependent, and the single-point formula requires a less 
straightforward treatment. 12 ' 13 As for the derivatives in 
Eq. (3), they can be evaluted in two ways: either "ana- 
lytically" (by means of perturbation theory), or "numer- 
ically" (by means of finite differences). 

The analytical-derivative approach relies on the per- 
turbation formula 

\d a u n0 ) = > \u m0 ) , (4) 

where v a is the a-component of the velocity operator 

v = i[H,r] = V k i/ k | k=0 . (5) 

Eq. (4) is convenient for tight- binding implementations, 
where the sum is over a small number of terms. We 
also notice that the matrix representation of r, for use in 
Eqs. (1) and (5), is usually taken to be diagonal on the 
tight-binding basis. 

The numerical-derivative approach looks more conve- 
nient for first- principle implementations, since it requires 
neither the (slowly convergent) sum over states, nor the 
matrix elements of the velocity operator. In order to give 
the most general formulation for non-rectangular cells it 
is expedient to switch to a coordinate-independent form. 
In order to shorten the equations, in all of the following 



developments of Eq. (3) a sum over the occupied states 
is implied. 

If hj are the shortest reciprocal vectors of the supercell, 
Eq. (3) can be identically recast as 

M ~ - £i3 ^y Im (d jUn0 \(H + e n0 2 M ) \d lUnQ ). 

(6) 

where a sum over ijl is implied, and dj indicates the 
partial derivative in the direction of hj : 

\djU n0 ) = lim — 1-| ( \u nX bj) - Ko) )• (7) 

Notice that Eq. (7) implicitly assumes that |u„k) is a dif- 
ferentiate function: but this is generally not the case 
when the eigenstates are obtained from numerical di- 
agonalization. The discretization must then be done 
in a specific gauge: as in I, we fix the problem by 
adopting the "covariant derivative" approach, introduced 
in Refs. 16 and 17. One defines the overlap matrix 
S'„„'(k) = (u„ |w„'k), and the "dual" states 

Kk) = ^5-, 1 „(k)K, k ) ) (8) 

n' 

which enjoy the property (u„o|itn'k) = &nn'- Using this, 
approximating Eq. (7) with its A = 1 value, and inserting 
in Eq. (6) we finally get 

M ~ - Im (u nbj |( H + 6» - 2 M ) \u nbl ). (9) 

Next, we wish to evaluate \u n bj) without actually diag- 
onalizing the Hamiltonian at k ^ 0. To this aim, we no- 
tice that the state e~ lhj ' r \u n0 ) obeys periodic boundary 
conditions and is an eigenstate of H n \ 3j corresponding, 
possibly, to a different occupied eigenvalue. The transfor- 
mation in Eq. (8) restores the correct ordering anyhow; 
therefore we can simply identify \u n bj) — c~* bj ' r |u„o), 
transform to the dual states by means of Eq. (8), and 
insert into Eq. (9) which eventually becomes the single- 
point, k = formula, aimed at. 

For a two-dimensional (2d) system the magnetization 
is a pseudoscalar, and the analogue of Eq. (9) reads 

M = J^IU^IU I Im («nbi|(-ffo + £«0 - 2/J, ) |u„b 2 ). 

c[Zn ) \ui | | U2| 

(10) 

Similarly, the single-point formula for the Chern number 
reads 

|bi x b 2 | 

= ~ 27r|bi| |b 2 | \ u nbi\unb 2 )- (11) 

As in previous works, 5 ' 6 we find expedient to val- 
idate the present findings on the Haldane model 
Hamiltonian: 18 it is comprised of a 2d honeycomb lat- 
tice with two tight-binding sites per primitive cell with 
site energies ±A, real first-neighbor hoppings t\, and 
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FIG. 1: Four unit cells of the Haldane model. Filled (open) 
circles denote sites with Eq = — A (+A). Solid lines connect- 
ing nearest neighbors indicate a real hopping amplitude t\; 
dashed arrows pointing to a second-neighbor site indicates a 
complex hopping amplitude t^e 1 ^ '■ Arrows indicate sign of the 
phase (j> for second-neighbor hopping. 



complex second-neighbor hoppings t%e f , as shown in 
Fig. 1. Within this two-band model, one deals with insu- 
lators by taking the lowest band as occupied. Following 
the original notations 18 we choose the parameters A = 1, 
t\ = 1 and |t 2 1 = 1/3. As a function of the flux parame- 
ter (j), this system undergoes a transition from zero Chern 
number to |C| = 1 when | sin 0| > 1/V3. Here we ad- 
dress periodic supercells made of L x L primitive cells, 
up to L = 32 (2048 sites), taking the lowest L 2 orbitals 
as occupied. 

Before actually addressing magnetization, we start 
benchmarking the accuracy of our single-point formula 
for the Chern number, whose value is known exactly as 
a function of the parameters of the model. The con- 
vergence of the Chern number — computed from Eq. (11) 
and its analytical-derivative analogue — as a function of 
the supercell size, is shown in Fig. 2, for <ft — OAtt, where 
the exact value is 1. Both approaches (analytical and nu- 
merical derivative) converge very fast. For instance the 
numerical-derivative approach yields an error of 7 x 10~ 3 
for L = 6, and smaller than 10~ 5 for L = 32. We are 
showing here the results for a <fi value well inside the 
C = 1 domain. We also find that the convergence wors- 
ens near the transition point | sin</>| = l/y/3. 

Numerical evaluation of Chern numbers is a staple tool 
in the theory of the quantum Hall effect, where super- 
cells are routinely used to account for disorder and/or 
electron-electron interaction. However, even in a super- 
cell framework, a discrete reciprocal mesh (or equiva- 
lently a mesh of phase boundary conditions) has been in- 
variably used in the algorithms implemented so far. 19-22 
Here we have shown that, provided the supercell is large 
enough, no mesh is needed: the Chern number can 
be evaluated from a single Hamiltonian diagonalization 
(with a single choice of boundary condition) . The ratio- 
nale behind our finding is simple: the Chern number is by 
definition an integral, whose integration domain shrinks 
to a single point in the limit of a large supercell. 

The single-point orbital magnetization M of the model 
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FIG. 2: Convergence of the Chern number as a function of 
the supercell size, evaluated using the single-point formulas 
(see text), for the Haldane model Hamiltonian at cj> = 0.4-7T. 
The largest L corresponds to 2048 sites. 



system, computed from Eqs. (3) and (10) as a function of 
the supercell size, is shown in Fig. 3, again for <p = OAtt. 
In this case the analytical-derivative approach converges 
definitely better, showing in fact the same kind of rel- 
ative error as the Chern number, while the numerical- 
derivative approach proves somewhat less accurate. 
In conclusion, we provide here the key formulas for 
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FIG. 3: Convergence of the orbital magnetization as a func- 
tion of the supercell size, evaluated using the single-point 
formulas (see text), for the Haldane model Hamiltonian at 
4> = 0.4-7T, The largest L corresponds to 2048 sites. 
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computing the orbital magnetization of a condensed sys- 
tem from first principles in a supcrccll framework and us- 
ing a single k point, to be used as they stand within Car- 
Parrinello simulations in an environment which breaks 
time-reversal symmetry. We have validated the present 
formulas on a simple tight-binding model Hamiltonian in 
2d, and checked their (fast) convergence with the super- 
cell size. Last but not least, we have proved that even 
the Chern number — which has a paramount relevance in 
quantum-Hall-effect simulations — can be computed from 
a single Hamiltonian diagonalization, and converges fast 
with the supcrccll size. 

We acknowledge fruitful discussions with D. Vanderbilt 
and T. Thonhauser. Work partly supported by ONR 
through grant N00014-03-1-0570. 

Appendix: More general boundary conditions 

The single-point formulas discussed so far are based 
on Eq. (7), with A = 1, and eventually require diagonal- 
izing the Hamiltonian at the T point only, ergo solving 
the Schrodingcr equation with periodic boundary con- 
ditions on the supercell. This is by far the most com- 
mon choice among Car-Parrinello practitioners, although 
other choices are possible. 

In order to extend our single-point formulas to more 
general boundary conditions it would be enough to switch 
from Eq. (7) (at A = 1) to alternative expressions for the 
directional derivatives. The only important requirement 
is that the two eigenstates therein differ by a supercell 



reciprocal vector. 

For the sake of simplicity we explicitly deal here only 
with the 2d case of antiperiodic boundary conditions, cor- 
responding to a zone-boundary single point: in fact, an- 
tiperiodic eigenstates obtain by choosing the special k- 
vector «i = (bi + b 2 )/2. It is then expedient to define 
even k 2 = (bi — b 2 )/2 = K\ — b 2 and to switch from 
Eq. (7) to 

|<9jU„o) = lim 2A ^..| ( K \Kj ) - \u n — AAtj ) ), (12) 

where now the j subscript indicates the derivative in the 
direction of Kj . In terms of such derivatives the magne- 
tization formula, Eq. (10) reads 

M = ~^Tv2\ — m — T Im (diu n0 \(H + e„ - 2^ ) \d 2 u n0 ), 
c(2tt) z \ki\ \k 2 \ 

(13) 

and similarly for the Chern number. 

In the case of a large supercell we approximate Eq. (12) 
with its A = 1 value, noticing that all the needed 
states obtain from a single Hamiltonian diagonalization 
at «x. In fact \u nK2 ) = c lh2 ' r \u nKl ) and \u n _ Kj ) = 
e i2Kj-r| Un ^ w \ iCIC 2kj is a reciprocal-lattice vector. 

While the states \u n «j) are to be used as they stand, 
the states \u n - Kl ), \u nK2 ), and \u n _ K2 ) must be reg- 
ularized to their dual counterpart, by means of the ob- 
vious analogue of Eq. (8). One then uses these states in 
Eq. (12) with A = 1 and finally in Eq. (13). 
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